Ejemplo de regresión completo

En este lab veremos qué es lo que tenemos que hacer para ajustar un modelo de regresión lineal múltiple y llegar a conclusiones finales

Datos: Accidentes de tráfico

Consideremos los datos Highway. Estos datos corresponden a un artículo de maestría no publicada de Carl Hoffstedt. Se relacionan a la tasa de accidentes de auto (en accidentes por millón de millas de vehiculos) y su asociación con varios términos potenciales. Los datos incluyen 39 secciones de Highways americanos grandes en el estado de Minnesoya en 1973. El objetivo original de estudio fue entender el impacto de variables de diseño: acpts, slim, Sig y shld que estaban bajo control del Departamento de Highways, en los accidentes.

Las variables en los datos son:

  • adt = conteo de tráfico diario promedio en miles
  • trks = volumen de camiones como porcentaje del volumen total
  • lane = número total de lineas de tráfico
  • acpt = número de puntos de acceso por milla
  • sigs = número de intercambios señalizados por milla
  • itg = número de tipos de intercambio de caminos por milla
  • slim = límite de velocidad en 1973
  • len = longitud de segmentos de Highway en millas
  • lwid = ancho de lines, en pies
  • shld = ancho del acotamiento en pies
  • htype = Una variable indicadora del tipo de camino o la fuente de fondeo del camino:
    • “mc” para colectores mayores,
    • “fai”para highways interestatales Federales
    • “pa” para highways de arteria principal
    • “ma” para highways arteriales mayores
  • rate tasa de accidentes de 1973 por millón de millas de vehículos
head(Highway)
##   adt trks lane acpt sigs  itg slim   len lwid shld htype rate
## 1  69    8    8  4.6 0.00 1.20   55  4.99   12   10   fai 4.58
## 2  73    8    4  4.4 0.00 1.43   60 16.11   12   10   fai 2.86
## 3  49   10    4  4.7 0.00 1.54   60  9.75   12   10   fai 3.02
## 4  61   13    6  3.8 0.00 0.94   65 10.65   12   10   fai 2.29
## 5  28   12    4  2.2 0.00 0.65   70 20.01   12   10   fai 1.61
## 6  30    6    4 24.8 1.84 0.34   55  5.97   12   10    pa 6.87
str(Highway)
## 'data.frame':    39 obs. of  12 variables:
##  $ adt  : int  69 73 49 61 28 30 46 25 43 23 ...
##  $ trks : int  8 8 10 13 12 6 8 9 12 7 ...
##  $ lane : int  8 4 4 6 4 4 4 4 4 4 ...
##  $ acpt : num  4.6 4.4 4.7 3.8 2.2 24.8 11 18.5 7.5 8.2 ...
##  $ sigs : num  0 0 0 0 0 1.84 0.7 0.38 1.39 1.21 ...
##  $ itg  : num  1.2 1.43 1.54 0.94 0.65 0.34 0.47 0.38 0.95 0.12 ...
##  $ slim : int  55 60 60 65 70 55 55 55 50 50 ...
##  $ len  : num  4.99 16.11 9.75 10.65 20.01 ...
##  $ lwid : int  12 12 12 12 12 12 12 12 12 12 ...
##  $ shld : int  10 10 10 10 10 10 8 10 4 5 ...
##  $ htype: Factor w/ 4 levels "mc","fai","pa",..: 2 2 2 2 2 3 3 3 3 3 ...
##  $ rate : num  4.58 2.86 3.02 2.29 1.61 6.87 3.85 6.12 3.29 5.88 ...

Paso 1: ve los datos

Lo primero que hay que hacer es examinar la gráfica de dispersión de puntos, y seleccionar algunas transformaciones a los predictores originales con la intención de (1) lograr normalidad de los datos

hw <- Highway %>% select(-htype)
pairs(hw,col=Highway$htype) # No se grafica la variable categórica "htype"; se usa para marcar los datos

summary(Highway)
##       adt             trks             lane            acpt      
##  Min.   : 1.00   Min.   : 6.000   Min.   :2.000   Min.   : 2.20  
##  1st Qu.: 5.00   1st Qu.: 8.000   1st Qu.:2.000   1st Qu.: 6.95  
##  Median :13.00   Median : 9.000   Median :2.000   Median :10.30  
##  Mean   :19.62   Mean   : 9.333   Mean   :3.128   Mean   :12.16  
##  3rd Qu.:24.00   3rd Qu.:11.000   3rd Qu.:4.000   3rd Qu.:14.60  
##  Max.   :73.00   Max.   :15.000   Max.   :8.000   Max.   :53.00  
##       sigs             itg              slim         len        
##  Min.   :0.0000   Min.   :0.0000   Min.   :40   Min.   : 2.960  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:50   1st Qu.: 7.995  
##  Median :0.0900   Median :0.1300   Median :55   Median :11.390  
##  Mean   :0.4005   Mean   :0.2964   Mean   :55   Mean   :12.884  
##  3rd Qu.:0.5800   3rd Qu.:0.3600   3rd Qu.:60   3rd Qu.:17.800  
##  Max.   :2.5100   Max.   :1.5400   Max.   :70   Max.   :40.090  
##       lwid            shld        htype         rate      
##  Min.   :10.00   Min.   : 1.000   mc : 2   Min.   :1.610  
##  1st Qu.:12.00   1st Qu.: 4.000   fai: 5   1st Qu.:2.630  
##  Median :12.00   Median : 8.000   pa :19   Median :3.050  
##  Mean   :11.95   Mean   : 6.872   ma :13   Mean   :3.933  
##  3rd Qu.:12.00   3rd Qu.:10.000            3rd Qu.:4.595  
##  Max.   :13.00   Max.   :10.000            Max.   :9.230

Algunas observaciones:

  1. La variable sigs tiene muchos 0s, que corresponden a los caminos que no tienen intercambios señalizados. Entonces se puede sustituir la variable transformando por \((sigs1 = sigs*len+1)/len\) para que tome valores positivos que pueden transformarse con potencias

  2. Podemos ver que las variables lwid y lane, aun sin ser categóricas, toman sólo algunos valores (la primera sólo valores en 11-13 y la segunda valores 2-4-6-8). Con poca variabilidad, pueden explicar muy poco de la respuesta. Podemos ya sea considerarlas factores, o removerlas del análisis

  3. Se mecionó en clase que las variables \(v\) que tienen \(max(v)/min(v) \geq 10\) se pueden transformar a logaritmos para escalar los números. En particular adt acpt y len cumplen la condición.

  4. Podemos ver en el último renglón de la matriz, que hay algunas relaciones marginales modestas con la respuesta.

  5. Muchos de los predictores tienen algún tipo de asociación entre ellos.

Para atender las observaciones anteriores, hacemos las siguientes transformaciones y volvemos a graficar

hw <- hw %>% mutate(sigs1 = (sigs*len+1)/len,
                    log_adt = log(adt),
                    log_acpt = log(acpt),
                    log_len = log(len)) %>%
             select(-sigs,-acpt,-len,-lwid,-lane,-adt)
hw <- hw[,c(5,1:4,6:9)]  #reordena las variables en el dataframe
pairs(hw,col=Highway$htype)

Paso 2: Aplica regresión “gráfica”

Les voy a dar las ideas básicas desde un punto de vista aplicado: Supongan que la verdadera relación entre la variable de respuesta y los predictores a través de la siguiente relación: \[ E(y|X) = g(\beta'X) \] para alguna función \(g\) (completamente desconocida y no especificada). Si esto es cierto y \(y\) depende de \(X\) sólo a través de la combinación lineal de \(y\) vs \(\beta'X\), entonces la gráfica de \(\{y,\beta'X \}\) nos puede dar un indicio de quién es \(g\). Las condiciones para que este resultado (Li y Duan 1989) funcione son las siguientes:

  1. Si los predictores están linealmente relacionados, o aproximadamente la relación entre los predictores es lineal en el scatterplot, o mejor aun, son mnormales, entonces podemos usar la regresión como la combinación lineal que necesitamos.

  2. Dado lo anterior, las transformaciones de los predictores deben buscar linealizar ó normalizar la relación entre predictores. Esta tarea puede ser laboriosa sin herramientas interactivas, porque hay que hacer muchas gráficas.

Por ejemplo, en el scatterplot de los datos que tenemos, particularmente las variable itg y sigs1 al relacionarse con las otras variables muestra patrones no lineales. Podemos intentar transformando para ver si mejora la relación lineal:

hw <- hw %>% mutate(itg2 = sqrt(itg),  #experimentando llegué a esta tarnsformación
                    log_sigs1 = log(sigs1),
                    htype = Highway$htype) %>%
            select(-itg,-sigs1)
pairs(hw, col=Highway$htype)

Paso 3: Ajusta el modelo de regresión y grafica \(y\) vs \(\hat{y}\).

En este paso podemos ver cuál sería la función \(g\) a escoger, lo que sugeriría la transformación de la variable de respesta. Si la gráfica muestra una clara tendencia no lineal, entonces la respuesta debe ser transformada para corresponder a esa relación no lineal.

m1 <- lm(rate ~ ., data=hw)  # incluye al factor que habíamos sacado. Estos factores usualmente no se transforman
summary(m1)
## 
## Call:
## lm(formula = rate ~ ., data = hw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.91012 -0.62102 -0.02195  0.69901  1.75345 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 10.307788   4.061093   2.538   0.0172 *
## trks        -0.062485   0.092180  -0.678   0.5036  
## slim        -0.062432   0.061986  -1.007   0.3228  
## shld        -0.006614   0.108084  -0.061   0.9517  
## log_adt     -0.603072   0.406737  -1.483   0.1497  
## log_acpt     1.064998   0.521566   2.042   0.0510 .
## log_len     -0.914522   0.366275  -2.497   0.0189 *
## itg2         0.767352   0.998767   0.768   0.4490  
## log_sigs1    0.714586   0.261965   2.728   0.0111 *
## htypefai     1.214930   1.589018   0.765   0.4512  
## htypepa     -0.856320   1.058392  -0.809   0.4255  
## htypema     -0.200021   0.862725  -0.232   0.8184  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.028 on 27 degrees of freedom
## Multiple R-squared:  0.8095, Adjusted R-squared:  0.7319 
## F-statistic: 10.43 on 11 and 27 DF,  p-value: 4.178e-07
plot(hw$rate,fitted(m1))
lw1 <- loess(hw$rate ~ fitted(m1),span = 0.15)
j <- order(hw$rate)
lines(hw$rate[j],lw1$fitted[j],col="red")

La transformación de la respuesta puede ser sugerida por la transformación de Box-Cox. Podemos obtener la potencia \(\lambda\) del siguiente modo:

a <-boxCox(m1)

a$x[which(a$y==max(a$y))]
## [1] -0.02020202

Entonces la transformación sugerida para la respuesta es \(\lambdahat \approx -0.02\) y esta transformacion incluye en su intervalo de confianza al logaritmo, así que podemos transformar la respuesta y volver a ajustar el modelo:

hw <- hw %>% mutate(log_rate =  log(rate),
                    htype = Highway$htype)

m2 <- lm(log_rate ~ . -rate, data=hw)
summary(m2)
## 
## Call:
## lm(formula = log_rate ~ . - rate, data = hw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48644 -0.10418  0.00738  0.12390  0.37630 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  3.424825   0.969455   3.533   0.0015 **
## trks        -0.021560   0.022005  -0.980   0.3359   
## slim        -0.019220   0.014797  -1.299   0.2050   
## shld        -0.007655   0.025802  -0.297   0.7690   
## log_adt     -0.163968   0.097095  -1.689   0.1028   
## log_acpt     0.184751   0.124507   1.484   0.1494   
## log_len     -0.228973   0.087436  -2.619   0.0143 * 
## itg2         0.115995   0.238423   0.487   0.6305   
## log_sigs1    0.172316   0.062536   2.755   0.0104 * 
## htypefai     0.306300   0.379327   0.807   0.4264   
## htypepa     -0.218802   0.252657  -0.866   0.3941   
## htypema     -0.148325   0.205948  -0.720   0.4776   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2455 on 27 degrees of freedom
## Multiple R-squared:  0.8002, Adjusted R-squared:  0.7189 
## F-statistic: 9.833 on 11 and 27 DF,  p-value: 7.58e-07
plot(hw$log_rate, fitted(m2))
lw2 <- loess(hw$log_rate ~ fitted(m2),span = 0.15)
j <- order(hw$log_rate)
lines(hw$log_rate[j],lw2$fitted[j],col="red")

Este modelo parece ser suficientemente adecuado para ajustar, aunque el coeficiente de determinación es menor que el del modelo m1. Previo a la última transformación, quizá sería adecuado analizar si no hay outliers y puntos influenciales que estén alterando el ajuste. Podemos analizar los residuales como siguiente paso, en el modelo m1.

residualPlots(m1)

##            Test stat Pr(>|Test stat|)  
## trks          1.6715          0.10663  
## slim          1.0171          0.31850  
## shld         -0.1145          0.90970  
## log_adt      -0.1793          0.85909  
## log_acpt      1.7525          0.09148 .
## log_len       1.1042          0.27961  
## itg2         -1.0653          0.29656  
## log_sigs1     1.1115          0.27652  
## htype                                  
## Tukey test    1.4371          0.15069  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
residualPlots(m2)

##            Test stat Pr(>|Test stat|)
## trks          0.4448           0.6601
## slim         -0.7055           0.4868
## shld         -0.6013           0.5528
## log_adt       1.1134           0.2757
## log_acpt      0.0399           0.9685
## log_len       0.3596           0.7220
## itg2         -0.1911           0.8499
## log_sigs1     0.3652           0.7179
## htype                                
## Tukey test   -0.8103           0.4178

Una variación de las gráficas de residuales, son las gráficas de modelos marginales. Esta función hace una prueba de especificación: agrega un término cuadratico a la relación entre los predictores y los residuales para verificar si la relación es curva o no.

residualPlots(m1)

##            Test stat Pr(>|Test stat|)  
## trks          1.6715          0.10663  
## slim          1.0171          0.31850  
## shld         -0.1145          0.90970  
## log_adt      -0.1793          0.85909  
## log_acpt      1.7525          0.09148 .
## log_len       1.1042          0.27961  
## itg2         -1.0653          0.29656  
## log_sigs1     1.1115          0.27652  
## htype                                  
## Tukey test    1.4371          0.15069  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
residualPlots(m2)

##            Test stat Pr(>|Test stat|)
## trks          0.4448           0.6601
## slim         -0.7055           0.4868
## shld         -0.6013           0.5528
## log_adt       1.1134           0.2757
## log_acpt      0.0399           0.9685
## log_len       0.3596           0.7220
## itg2         -0.1911           0.8499
## log_sigs1     0.3652           0.7179
## htype                                
## Tukey test   -0.8103           0.4178

También podemos qq-plots, que están adaptados al caso de RLM

qqPlot(m1,id.n=3) #identifica las tres observaciones con los residuales más grandes

## [1]  5 26
qqPlot(m2,id.n=3)

## [1] 34 36

Por último, podemos anañizar los posibles puntos infuenciales:

influenceIndexPlot(m2)

influencePlot(m2) #combina los residuales estudentizados con distancias de Cook y valores de la matriz sombrero.

##       StudRes       Hat      CookD
## 25 -1.7252819 0.3423982 0.12034379
## 34 -2.5622986 0.2785552 0.17514351
## 36 -1.9192166 0.1934702 0.06697469
## 38  0.7499658 0.5857165 0.06735770
## 39 -0.7499658 0.5857165 0.06735770

Al parecer podemos intentar rehacer el análisis eliminando la observación 34, y posiblemente tengamos un mejor ajuste

hw34 <- hw[-34,]
m2a <- lm(log_rate ~ . - rate , data=hw34)
summary(m2a)
## 
## Call:
## lm(formula = log_rate ~ . - rate, data = hw34)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42397 -0.08137 -0.01081  0.13773  0.37102 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.546381   0.884011   4.012 0.000453 ***
## trks        -0.006385   0.020893  -0.306 0.762333    
## slim        -0.025276   0.013679  -1.848 0.076048 .  
## shld         0.016735   0.025349   0.660 0.514947    
## log_adt     -0.167700   0.088422  -1.897 0.069046 .  
## log_acpt     0.153026   0.114044   1.342 0.191259    
## log_len     -0.198476   0.080500  -2.466 0.020598 *  
## itg2         0.101466   0.217171   0.467 0.644237    
## log_sigs1    0.200212   0.057973   3.454 0.001908 ** 
## htypefai     0.227898   0.346749   0.657 0.516802    
## htypepa     -0.300035   0.232231  -1.292 0.207734    
## htypema     -0.075962   0.189640  -0.401 0.692018    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2235 on 26 degrees of freedom
## Multiple R-squared:  0.8315, Adjusted R-squared:  0.7602 
## F-statistic: 11.67 on 11 and 26 DF,  p-value: 1.886e-07
plot(hw34$log_rate, fitted(m2a))

boxCox(m2a)

En este último modelo podemos ver que algunos predictores no son significativos, por lo que podemos evaluar eliminarlos, considerando el ajuste de todos los modelos con estos predictores para obtener el mejor.

step(m2a)
## Start:  AIC=-104.29
## log_rate ~ (rate + trks + slim + shld + log_adt + log_acpt + 
##     log_len + itg2 + log_sigs1 + htype) - rate
## 
##             Df Sum of Sq    RSS      AIC
## - trks       1   0.00467 1.3036 -106.154
## - itg2       1   0.01091 1.3098 -105.973
## - shld       1   0.02177 1.3207 -105.659
## <none>                   1.2989 -104.291
## - log_acpt   1   0.08995 1.3888 -103.746
## - slim       1   0.17056 1.4695 -101.602
## - log_adt    1   0.17970 1.4786 -101.367
## - log_len    1   0.30369 1.6026  -98.307
## - htype      3   0.52532 1.8242  -97.384
## - log_sigs1  1   0.59584 1.8947  -91.943
## 
## Step:  AIC=-106.15
## log_rate ~ slim + shld + log_adt + log_acpt + log_len + itg2 + 
##     log_sigs1 + htype
## 
##             Df Sum of Sq    RSS      AIC
## - itg2       1   0.01444 1.3180 -107.736
## - shld       1   0.03091 1.3345 -107.264
## <none>                   1.3036 -106.154
## - log_acpt   1   0.08757 1.3911 -105.684
## - slim       1   0.18064 1.4842 -103.223
## - log_adt    1   0.18514 1.4887 -103.108
## - log_len    1   0.32341 1.6270  -99.733
## - htype      3   0.54355 1.8471  -98.910
## - log_sigs1  1   0.68196 1.9855  -92.165
## 
## Step:  AIC=-107.74
## log_rate ~ slim + shld + log_adt + log_acpt + log_len + log_sigs1 + 
##     htype
## 
##             Df Sum of Sq    RSS      AIC
## - shld       1   0.02881 1.3468 -108.914
## <none>                   1.3180 -107.736
## - log_acpt   1   0.07856 1.3966 -107.536
## - log_adt    1   0.17997 1.4980 -104.872
## - slim       1   0.19074 1.5087 -104.600
## - log_len    1   0.35547 1.6735 -100.662
## - htype      3   0.84669 2.1647  -94.882
## - log_sigs1  1   0.72374 2.0417  -93.104
## 
## Step:  AIC=-108.91
## log_rate ~ slim + log_adt + log_acpt + log_len + log_sigs1 + 
##     htype
## 
##             Df Sum of Sq    RSS      AIC
## <none>                   1.3468 -108.914
## - log_acpt   1   0.12455 1.4714 -107.553
## - log_adt    1   0.16303 1.5098 -106.572
## - slim       1   0.18728 1.5341 -105.967
## - log_len    1   0.42026 1.7671 -100.594
## - htype      3   0.83505 2.1819  -96.582
## - log_sigs1  1   0.70358 2.0504  -94.943
## 
## Call:
## lm(formula = log_rate ~ slim + log_adt + log_acpt + log_len + 
##     log_sigs1 + htype, data = hw34)
## 
## Coefficients:
## (Intercept)         slim      log_adt     log_acpt      log_len  
##     3.27719     -0.01962     -0.13859      0.16805     -0.22216  
##   log_sigs1     htypefai      htypepa      htypema  
##     0.20563      0.29093     -0.28367     -0.09170

El modelo sugerido es el que tiene el AIC más pequeño.